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consequences  for  the  interpretation  of  the  behavior  of  copper- 
gold  alloys  are  indicated. 
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Abstract : 

We  describe  the  results  of  Monte  Carlo  computations  of  the 
coherent  phase  diagram  (in  the  temperature-composition  plane)  of 
ordering  binary  alloys  on  a  face-centered  cubic  lattice.  Results 
on  long-  and  short-range  order  parameters  as  well  as  ordering 
energies  are  also  given.  We  consider  the  system  with  nearest 
neighbor  interaction  in  the  grand-canonical  ensemble  (equivalent 
to  an  Ising  antiferromagnet  in  a  magnetic  field)  as  well  as  in 
the  canonical  ensemble  (fixed  composition).  The  close  agreement 
between  both  approaches  serves  as  a  severe  test  of  the  good  accuracy 
'•-.bteined,  and  so  does  the  fact  that  differences  between  results  for 
lattices  with  2048  sites  and  with  16384  sites  are  found  to  be  negli- 
niblv  small.  Results  with  next-nearest  nej.ohbor  interaction  are  also 
given,  and  for  both  models  a  comparison  with  other  available 
predictions  is  made,  particularly  with  the  cluster-variation 
method.  While  the  latter  is  found  to  be  quite  accurate  at  stoi- 
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chiometric  composition,  it  appears  to  do  less  well  in  the 
more  general  case.  The  tetrahedral  approximation  of  the 
cluster  variation  method  predicts  a  topology  of  the  phase 
diagram,  in  the  case  of  nearest-neighbor  interaction  different 
from  the  computer  simulations.  Some  consequences  for  the 
interpretation  of  the  behavior  of  copper-gcld  alloys  are 
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1  .  Introduction 

We  have  investigated  the  behavior  of  ordering  binary  alloys 

such  as  the  copper-gold  system  on  a  rigid  face-centered  cubic 

(fee)  lattice,  within  the  framework  of  the  Ising-model  for  binary 

alloys  (AB)  where  one  assumes  pairwise  interactions  between  the 

atoms.  Let  us  denote  the  local  concentration  at  lattice  site 

i  by  :  c^=0  if  it  is  occupied  by  an  A-atom,  c^=1  if  it  is 

occupied  by  a  B-atom,  and  denote  the  interaction  energy  of 

AA  BB  AB 

AA-,  BB-  and  AB  pairs  as  v^,  v^  and  ^  ,  respectively.  Dis¬ 
regarding  vacancies  or  any  other  lattice  defects  the  energy 
of  the  system  is  then  given  by 


v  ■  I  {( l-c. )( 1-c. )v^+[c.  ( I-c O+c.C  1-c .)  Jv^+c . c .  v?1? }  , 
ilfj  i  J  ij  i  J  j  i  ij  i  J  ij 


AB. 


BB! 


(1) 


the  sums  in  Eq.  (1)  are  taken  over  all  pairs  of  neighbors  once. 

As  is  well  known,  for  the  ordering  behavior  of  this  system  it 

AA  BB  AB 

is  only  the  ordering  energy  (vi  ^+v^  j~2vi  j )  which  matters,  and 

AA  BB  AB 

not  the  three  energies  v^,  vi ^  ,  v^  separately  /I/.  This  fact 

is  utilized  by  transforming  Eq .  (1)  to  the  equivalent  problem 

of  an  Ising  magnet  in  a  magnetic  field:  a  "spin  ud"  ("spin  down") 

at  a  site  i,  S^=+1(-1)  can  be  identified  with  an  A(B)  atom  at  that 

site,  by  using  S^=1-2c^.  The  Ising  hamiltonian  is  then  expressed 

AA  BB  AB 

in  terms  of  a  spin-spin  interaction  j=” j+vi j“2vi j ) /4  and 

a  magnetic  field  H  =  -1  E  (v^-v®®)  +  X  where  j  and  u 

2  j(4i)  i3  iD  2  A  B  A  B 

are  the  chemical  potentials  of  the  A  and  B  atoms: 


*  m  -  ttj  ‘ H  ' 

Clearly,  for  general  ranae  of  the  interaction  the  description 

of  the  ordering  behavior  of  the  model  described  by  Eq .  (2)  is  a 
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formidable  problem.  Hence  we  restrict  the  range  of  the  inter¬ 
action  to  nearest  and  next  nearest  neighbors  on  the  lattice, 

as  is  usually  assumed  /2-4/.  Denoting  nearest  neighbor  and  next- 
nearest  neighbor  pairs  by  <nn>  and  <nnn->,  respectively,  we  have 

K  =  -J  I  S  S.  +  aj  r  S  S,  -HIS.,  (3) 

<nn>  3  <nnn>  1  -1  i 

where  the  parameter  -a  measures  the  relative  strength  of  the  next- 
nearest  neighbor  interaction.  This  choice  of  a  model  is  also 
reasonable  if  one  tries  to  fit  real  alloys,  like  the  copper- 
gold  system,  in  terms  of  a  pairwise  interaction  model,  Eq . 

(1)  /5,6/.  There  is,  however,  at  present  no  general  agreement  on 
whether  this  model  (with  parameters  J,a  smoothly  depending  on 
the  relative  composition  cB  of  the  alloy,  Cg  being  the  thermal 
average  <ci>T)  is  appropriate,  or  whether  one  needs  a 
model  with  interactions  between  nearest  neighbor  pairs,  triplets 
or  quadruplets  of  atoms  /7,8/.  It  might  even  be  necessary  to  use 
explicitly  the  long-range  interaction  due  to  the  conduction 
electrons  of  the  metallic  alloy  /9/.  In  our  opinion  a  consid¬ 
erable  part  of  this  uncertainty  about  an  appropriate 
model  for  the  ordering  behavior  of  alloys  is  due  to  the  fact 
that  the  properties  of  these  models  are  not  known  with  the 
precision  needed  to  distinguish  among  them.  As  a  step 
in  this  direction,  the  present  paper  oresents  a  study  of 
the  model,  Eq.  (3)  ,  making  use  of  the  well-known  Monte  Carlo 
computer  simulation  method  / 10/.  The  phase  diagram  of  this  model 
system  in  the  T-a  plane  at  equal  concentrations  of  A  and  B 
atoms  (corresponding  to  H*0)  has  already  been  reported  /II/. 

Some  results  at  general  composition  have  also  been  presented 
very  briefly  /1 2/  for  -he  case  a»0  with  computations  done  in 
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the  grand-canonical  ensemble  (H  fixed  as  an  independent  variable) 
Here  we  present  more  extensive  results  of  this  latter  work, 
supplementing  it  with  computations  done  at  other  lattice  sizes, 
computations  in  the  canonical  ensemble  (c^  fixed  as  independent 
variable)  ,  and  preliminary  results  for  ci4=0.  It  turns  out  that 
even  the  nearest  neighbor  case  (ot  =  0,  H  =f  0)  is  of  great 
complexity:  thanks  to  a  large  body  of  work  /1 2-2  7/  the 
ground  states  of  this  model  are  quite  well  understood,  but  for 
nonzero  temperature  the  approximation  schemes 
available  so  far  give  greatly  varying  results  /1 8-23/. 

As  we  shall  see  below,  none  of  these  schemes  is  completely 
satisfactory. 

In  Section  2  we  give  a  description  of  the  different 
computational  procedures  for  the  quantities  of  interest  (order 
parameters,  ordering  energies, phase  diagram).  In  Section  3 
our  results  for  a=0  are  presented  and  compared  to  other 
predictions  /17-23/.  Section  4  presents  our  results  for  a=-0.25 
and  compares  them  to  the  cluster  variatio  (CV)  predictions  /4/,' 
while  Section  5  summarizes  our  conclusions. 


2.  The  Monte  Carlo  Method 


We  consider  a  system  of  N=4L^  spins  on  an  fee  lattice 
with  periodic  boundary  conditions.  For  the  study  of  ordering 
phenomena,  L  has  to  be  chosen  such  that  the  periodic  boundary 
condition  is  compatible  with  the  ordered  structures  expected. 

Fig.  1  shows  some  of  the  orderings  which  have  been  identified 
as  ground-state  orderings  for  our  model  /2,3/,  translating 
them  into  the  spin  representation.  While  the  AB  structure  of 
CuAuI-type  as  well  as  the  A^B  structure  of  Cu^Au-type  (which 
both  occur  for  a  ^  0)  fit  the  periodic  boundary  condition 
for  any  L,  both  the  A2B2  structure  and  the  A3B  structure  of  Al^Ti 
type  (which  occur  for  a<0)  fit  the  periodic  boundary  condition 
only  for  even  L.  For  odd  L  this  structure  would  necessarily  require 
an  interface  ("domain  wall"),  which  would  disturb  the  results.  Hence 
L=8  and  L=16  ( N— 1 6384)  were  chosen  to  study  these  structures,  and 
to  make  sure  that  finite  size  effects  are  already  negligible 
(early  work  /24/  using  L=5,  i.e.  N=500  and  very  short  "observation 
times"  could  not  study  reliably  the  phase  transitions,  due  to 
f:nite  size-  and  finite  time-effects) .  finally,  for  the  A2B- 
structure,  L  has  to  be  a  multiple  of  3  in  order  to  avoid  domain 
effects,  as  is  obvious  from  Fig.  1e,  and  hence  L=15  (i.e.,  N=13500) 
was  chosen  for  studying  that  structure. 

Similar  considerations  are  possible  for  other  structures 
which  occur  in  the  ground  state  /2,3/,  such  as  the  AjB,A4B 
and  AgBj  structures.  However,  no  attempt  was  made  to  study 
these  structures  at  finite  temperatures  in  the  simulations: 
in  the  temperature  recrime  of  interest  here  (k_T/|j|  >  0.3) 
even  the  A,B  structure  was  never  found  tc  be  stable,  and  we 
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expect  a  similar  situation  for  these  other  more  complicated 
structures,  as  will  be  discussed  below. 

In  the  case  of  a^O  the  occurrence  of  long  period  super¬ 
structures  (such  as  the  CuAuII  structure)  is  possible,  and 
in  other  problems  of  Ising  models  with  competing  interactions 
such  incommensurate  structures  have  in  fact  been  identified 
/25-27/.  Since  the  periodic  boundary  condition  is  inconsistent 
with  such  structures  for  any  finite  L,  the  analysis  of  this 
case  by  Monte  Carlo  methods  is  rather  delicate  /25/.  We  do 
not  consider  this  possibility  here,  since  both  in 
the  theoretical  models  /25-27/  and  in  the  experiments  on  the 
copper-gold  system  /28/  the  regime  of  stability  of  the  in¬ 
commensurate  phase  in  the  phase  diagram  is  very  small. 

The  microscopic  state  of  the  system  is  specified  by  the 
configuration  X, X={S1 ,S2 , . . . ,SN) .  When  the  system  is  at 
equilibrium,  the  probability  of  a  configuration  X  is 

Peq(X)  exp[-H,(X)/kBT]  ,  Z  =  l  exp[-«  (X)/kBT]  ,  (4) 

where  is  Boltzmann^  constant  and  T  the  temperature  of  the 
system,  2£(X)  being  given  by  Eq.  (3).  Starting  from  an  initial 
configuration  XQ  (the  choice  of  which  is  discussed  below) ,  one 
lets  the  configuration  of  the  system  evolve  according  to  the 
following  algorithm  /10/:  using  pseudorandom  numbers  one  generates 
a  random  change  of  the  configuration  $-»X'.  In  the  case  of  the 
grand-canonical  ensemble  {  magnetic  field  in  Eq.  (3)  is  held 
fixed,  concentration  c0  is  allowed  to  fluctuate}  this  transition 
X-X'  is  taken  to  be  the  flip  of  a  randomly  chosen  spin. 
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In  the  case  of  the  canonical  ensemble,  a  nearest  neighbor 
AB-pair  is  chosen  at  random  and  interchanged.  This  "Kawasaki 
dynamics”  /29/  conserves  the  concentration  Cg  of  the  system. 

As  is  well  known,  both  ensembles  should  yield  results  equivalent 
to  each  other  in  the  thermodynamic  limit  N-®.  Since  we  work 
with  finite  N  it  is  important  to  see  that  the 
results  do  not  depend  on  the  ensemble  chosen,  at  least  within 
reasonable  error  limits. 

In  both  cases  one  computes  the  energy  change  6u«  '&U%' )  -  ^t(X) 
resulting  from  the  configurational  change.  The  "transition  probability"  / 10/ 
W=exp(-6u/k  _T)/[1+exp(-6u/knT) ]  is  then  compared  to  a  random 

D  D 

number  n»  chosen  uniformly  from  the  interval  [0,1].  If  W  >  n 
the  transition  is  performed;  otherwise  the  old  configuration 
X  is  counted  once  more  for  the  averaging,  the  attempted  X'  is 
rejected, and  another  transition  is  tried. 

This  algorithm  is  known  /10/  to  lead  to  the  thermal 
equilibrium  distribution,  Eq.  (4),  in  the  limit  where  the 
number  of  states  generated  tends  to  infinity.  In  practice,  one 
hopes  to  achieve  accurate  results  when  this  number  lies  in  the  range  from  a  few 
hundred  to  a  few  thousand  Monte  Carlo  steps/site  (MCS) .  To  accomplish 
this  aim,  it  is  important  to  choose  the  initial  states  XQ 
appropriately.  For  runs  using  the  grand-canonical  ensemble, 
it  is  most  convenient  to  use  the  fully  ordered  states  (Fig.  1) 
as  the  initial  condition.  Although  the  ground-state  is  known 
to  be  highly  degenerate  for  a«0  /12-17/,  it  is  known  that  at 
T>0  only  the  most  symmetric  structures  contribute  /30/.  Thus 
for  a«0  only  the  structures  of  Fig.  1a,c,e  were  used,  while 
for  a*-0.25  we  used  the  structures  of  ?io.  1b, c.  As 


we  shall  see  below,  the  transition  from  the  ordered  phases  to 
the  disordered  one  are  of  first  order.  Consequently,  one  r.as 
to  pay  attention  to  the  occurrence  of  hysteresis:  the  limit 
of  stability  of  the  various  ordered  phases  is  not  identical 
with  the  location  of  the  first-order  phase  transition,  but  rather 
locates  the  regime  where  the  lifetime  of  metastable  states  becomes 
short  due  to  nucleation  processes.  Hence  we  also  used  states 
with  a  random  or  a  ferromagnetic  spin  configuration  as  appropriate 
initial  states  for  the  disordered  phase.  As  it  is  more 
difficult  to  create  order  out  of  disorder  than  vice  versa,  ne 
finds  even  more  pronounced  metastability  effects  when  one 
approaches  the  first-order  transition  from  the  disordered 
But  by  a  careful  analysis  of  the  relaxation  of  the  metasta. -e 
states  on  both  sides  of  the  transition  one  can  locate  the  phase 
boundary  with  reasonable  accuracy  / 3 1 / -  This  will  become  apparent 
from  the  detailed  examples  shown  in  Section  3. 

For  runs  using  the  canonical  ensemble,  where  the  concentra¬ 
tion  c  is  kept  at  that  of  the  initial  state,  the  fully  ordered 
configurations  can  be  used  as  initial  states  only  for  the 
respective  stoichiometric  compositions.  Non-stoichiometric 
initially  ordered  states  could  be  produced  by  randomly  changing 
the  appropriate  fraction  of  A-atoms  into  B-atoms  (or  vice  versa) 
in  such  a  fully  ordered  state.  Most  of  our  runs  in  the 
canonical  ensemble  were  started  from  a  random  initial  configuration. 

We  conclude  this  section  by  summarizing  the  definitions  of 
the  order  parameters  which  we  compute,  paying  attention  to  the 
appropriate  translation  from  one  ensemble  to  the  other.  Labelling 
the  four  sublattices  of  the  unit  cell  of  the  fee  lattice  as 


10 


indicated  in  Fig.  1,  we  introduce  the  sublattice  magnetizations 

m  , 

v 


m  =  (1/N)  I  <S .  >  ,  v=1  ,2  ,3  ,4  (5) 

i€  v  1 

where  <  >  is  computed  by  taking  time  averages  over  one  or  more 

Monte  Carlo  runs  /lo/.  Then  the  order  parameter  components  of  the 
AB-structure  (Fig.  la)  are  expressed  in  terms  of  the  as 


~( 1 ) 
raA8 


~(2)  ~(3) 

®AB  ml  m2  ®3+m4 ’  mAB 


cn | 


(6) 


while  the  order  parameter  components  of  the  A^B-structure  (Fig.  1c) 
are 


~(l)  ~(2)  ~(3)  ~(4) 

mA3B*ml  +m2+m3-m4  *  mA3B”Vm2_m3+:V  mA3B’,VVIVm4  *  mA3B“‘ml  *Wm4  (7) 


Noting  that  c^U-S^/2,  the  average  magnetization  m. 


m=m1 +m2+m3+m4  , 


(8] 


is  related  to  the  average  concentration  c0  of  B  atoms  in  the  alloy 
by 


c0  =  (1-m)/2  .  (9) 

The  order  parameter  m  {Eqs.  (6), (7)}  is  related  to  the  standard 
long-range  order  parameter  /I/  &  by 

’()  =  (m-m)/(1-m)  .  (10) 

Short  range-order  in  the  crand-canonical  ensemble  is  described 


1 1 


in  terms  of  the  correlation  functions  [R  being  the  distance 
between  sites  i,j  on  the  lattice] 

g (R)  =  <S  S  ->  -  m2  (11) 

This  quantity  is  related  to  the  standard  Cowley  short  range  order 
parameters  a(R)  /I ,32/  by 

a<3)  =  g ( R) /  ( 1 -m2  >  (12) 

In  this  paper,  we  only  consider  for  R  a  nearest- 

neighbor  distance. 

Finally  we  consider  the  energy  u,  which  in  the  grand- 
canonical  ensemble  is  defined  as  follows 


u  =  <  Jt  >/N  ■  U  .  -mH  ; 

int 


(13) 


this  quantity  is  related  to  the  usual  ordering  energy  /2/ 

AU  =  I  i  (2v^-v^-v^ )  [  c .  ( 1  -  c  . )  +  c  •  (  1  -  c  .  )  ]/N  as 
^  j  ^  1  J  1  J  1  J  1  J  J  1 


AU  -  fUint  +  J<°) J  ' 


(14) 


where  J (O)  =  Z  J../N.  We  note  from  Eqs.  (11)— (14)  that  the 

i+3  3 

ordering  energy  can  be  expressed  in  terms  of  the  Cowley  short 
range  order  parameters  as 


AU  =  4cQ ( 1 -c_) ( J (0) -  Z  J ,  .  a  (  R)  /N  ]  .  (15) 

D  D  .  ,  .  1  j 
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3 .  Results  for  the  Mearest-Neicihbor  Interaction  fq=0) 

Figs.  2-4  show  typical  examples  of  "raw  data"  of  the  Monte 
Carlo  computation.  Most  of  the  data  were  obtained  in  the  grand- 
canonical  (GC)  ensemble,  N=16384.  It  turned  out  that  away  from 
the  transitions  it  was  sufficient  to  average  over  180  MCS/site, 
after  emitting  the  initial  60  MCS/site  of  the  simulation  to  get 
rid  of  the  influence  of  the  initial  condition,  while  close  to 
the  transition  averages  were  taken  over  900  MCS/site.  Runs 
using  the  canonical  (C)  ensemble  were  performed  at  the  compositions 
c  =0.225,  0.250,  0.265,  0.275,  0.30,  0.35,  0.40,  0.475  and  0.50, 

O 

respectively.  Identifying  the  corresponding  field  from  the  grand- 
canonical  m  vs.  H-curve  and  Eq.  (9),  the  energies  and  order 
parameters  of  the  C  and  GC  ensemble  could  be  compared.  Figs.  2-4 
show  that  very  good  agreement  between  these  two  types  of  calcula¬ 
tion  was  in  fact  obtained.  In  addition,  data  taken  for  different 
N  agree  with  each  other  well. 

From  Figs.  2-4  we  note  that  the  magnetization  varies  linearly 

with  the  field  in  the  disordered  phase  (  and  therefore  so  does  the  order 

parameter  m.  =m/2  in  this  range) .  In  the  regime  with  long-range 

order,  however,  the  magnetization  deviates  from  this  simple  linear 

law,  and  at  the  first-order  transitions  (occurring  at  one  of  the 

three  critical  fields  or  H^,  respectively)  one  can 

clearly  identify  the  jump  of  the  magnetization,  Am»m+-m_, 

where  m  (m  )=lim  m(H). Since  there  is  no  one-phase  equilibrium 
H-H  (H j 

state  with  m  <m<m. ,  it  is  clear  that  states  in  this  interval  are 

••  *r 

two-phase  mixtures,  and  hence  one  obtains  the  boundaries  of  the 
two-phase  regions  in  the  T-c  phase  diagram  from  using  m_,m+  in 
Eq.  (9) .  The  critical  fields  are  estimated  by  performing  runs 
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with  different  initial  conditions  and  searching  for  relaxation 
of  metastable  states,  as  mentioned  in  the  previous  section. 

Figs.  2,3  include  a  few  examples  of  points  where  the  A^B  ordering 
was  found  initially  metastable  while  longer  runs  revealed  a 
clear-cut  relaxation  towards  the  disordered  phase. 

Data  as  displayed  in  Figs.  2-4  have  been  taken  at  kgT/lJl® 

0.3,  0.4,  0.6,  0.8,  1.0,  1.2,  1.4,  1.5,  1.6,  1.7,  1.75,  1.8,  1.85, 
1.9,  and  2.0,  respectively,  and  thus  the  phase  boundaries  are 
constructed  both  in  the  T-c_  representation  (Fig.  5a)  and  in 

S 

the  T-H  representation  (Fig.  5b) .  We  think  that  the  relative 
accuracy  of  our  estimates  for  the  critical  fields  (which 

merge  at  Hcl  for  T-»0)  is  of  the  order  of  a  few  percent,  i.e.  the 
size  of  the  circles  in  Fig.  5.  The  accuracy  of  the  phase  boundaries 
in  the  T-c0  plane  is  better  by  about  a  factor  of  2,  because  the 
"susceptibility"  x  =  3nt/  OH/ [  J \  )  is  small.  Note  also  (Fig.  2), 
that  a  small  error  in  the  location  of  the  critical  fields  does 
not  affect  the  magnitude  of  the  jump  of  the  magnetization  very 
much,  since  (close  to  the  respective  critical  fields)  the  m  vs .  h- 
curves  are  nearly  parallel  in  the  disordered  phase  and  the  A^B  phase 
Therefore  we  think  that  meaningful  estimates  for  the  width  of  the 
two-phase  regions  have  been  obtained. 

Fig.  6  shows  the  temperature  variation  of  various  quantities 
at  fixed  field.  The  estimated  transition  temperatures  (taken 
from  Fig.  5)  are  included.  It  is  seen  that  pronounced  super¬ 
heating  of  the  ordered  phase  (and  supercooling  of  the  disordered 
phase)  occurs,  when  the  transition  temperature  is  low.  This  meta¬ 
stability  is  expected  from  the  shape  of  the  phase  boundaries  in 
Fig.  5,  of  course,  since  varying  the  temperature  at  either  fixed 
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field  or  fixed  composition  one  cuts  the  phase  boundary  typically 

at  rather  small  angles.  This  is  also  the  reason  why  in  this  case 

the  computations  in  the  canonical  ensemble  are  convenient  for 

locating  the  transition  only  for  compositions  near  stoichiometry. 

Due  to  these  hysteresis  phenomena,  no  attempt  v.as  iv.aJ< 

to  estimate  the  width  of  the  two-phase  region  from  the  runs  in 

the  canonical  ensemble.  But  apart  from  this  latter  ambiguity, 

both  calculations  give  the  same  phase  boundary,  as  is  seen  in 

Fig.  7  where  we  compare  our  results  to  the  phase  diagram 

predicted  by  the  cluster  variation  (CV)  method  /20,22/ .  It  is 

seen  that  the  CV  method  overestimates  the  transition  temperatures 

in  the  stoichiometric  cases  only  by  a  few  percent,  while  it  is 

very  inaccurate  for  intermediate  compositions.  In  addition,  it 

seems  to  underestimate  somewhat  the  widths  of  the  two-phase 

regions.  An  interesting  feature,  however,  which  is  correctly 

predicted  by  the  CV  method  is  that  the  maximum  transition 

temperature  of  the  A^B  phase  does  not  occur  at  the  stoichiometric 

composition  c0=1/4  but  rather  for  cQ  *  0.265.  It  should  also  be 

noted,  that  simpler  approximations  (like  the  Bragg-Wi lliams 

approximation,  Bethe-  and  quasichemical  approximations  /18,19,21/J 

fail  even  more  dramatically  than  the  CV  method,  as  they  predict 

not  only  a  wrong  topology  of  the  phase  diagram,  but  cannot  even 

locate  the  transition  temperatures  in  the  stoichiometric  cases 

with  reasonable  accuracy  (see  / 1 2/ ) .  On  the  other  hand,  the  method 

of  Wu  and  Tahir-Kheli  /6/  (which  refers  to  a  model  including 
nd  rd 

2  and  3r  nearest  neighbor  interaction)  predicts  all  transitions 
to  the  disordered  phase  to  be  of  second  order;  even  poorer  are 
the  real-space  renormalization  results  /23/  which  fail  to  yield 
the  ordering  of  the  AB  phase,  and  predict  the  transition  of  the 
A,B  phase  tc  be  c?  second  order.  As  »  result,  no  analytic  method 
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is  presently  capable  to  predict  a  satisfactory  phase  diagram 
for  this  model,  and  an  extension  of  the  methods  of  Refs.  6, 

22  and  23  seems  to  be  called  for. 

While  the  CV  method  predicts  a  triple  point  at  nonzero 
temperature  Tt<  below  which  one  can  observe  transitions  from 
the  A^B  structure  to  the  AB  structure,  we  have  tentatively 
extrapolated  our  results  (Fig.  5)  such  that  Tt=0.  Some 
justification  of  this  extrapolation  is  found  when  one  considers 
the  ground-state  energy  as  function  of  composition  (Fig.  8a) 
or  field  (Fig.  8b) .  It  turns  out  that  the  variation  with 
composition  is  much  more  complicated  to  understand  than  the 
behavior  as  function  of  the  field.  While  for  a  -*  0+  the  ground 
state  is  a  two  phase  mixture  of  pure  A  and  A^B  phases  for 
0  <  CB  <  ?  '  311(3  of  A3B»AB  phases  for  "J  <  cg  <  'J  »  several 
other  phases  occur  for  a  -*  0~  /2,3/.  One  first  has  a  transition  to  a 
(nonstoichiometric)  A^B  phase;  but  the  location  of  this  transition  is 
uncertain.  At  c0=1/6  one  starts  to  have  a  three-phase  mixture, 
containing  the  A^B , A^B  and  A^B  phases.  In  view  of  this  situation, 
it  is  interesting  to  note  that  the  outer  boundary  of  the  two- 
phase  region  between  the  disordered  and  the  A^B  phase  seems  to 
extrapolate  towards  cB=1/6(Fig.  5).  For  <  cQ  <  a  two- 
phase  mixture  between  A^B  and  AjB  is  predicted,  while  for  cB>1/3 
a  two-phase  mixture  between  AjB  and  A^B^  should  occur  which  has 
a  transition  (at  a  so  far  unknown  concentration)  to  a  pure  non¬ 
stoichiometric  AjB^  phase. 

Unfortunately,  we  have  not  been  able  to  establish  the 
connection  between  these  ground-state  phases  and  the  behavior 
seen  at  nonzerc  temperature .  In  par*  cul  sr ,  starting  the 
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system  in  the  A^B  structure  revealed  that  this  structure  was 
not  stable  (and  not  even  metastable)  down  to  k  T/  J  w  0.3, 

D 

and  we  expect  that  it  will  not  be  stable  at  any  nonzero 
temperature . 

The  reason  for  the  different-  instability  of  the  various 
ordered  phases  is  seen  more  clearly  in  Fig.  8b,  which  shows  that 
the  AB  and  A^B  structures  are  ground  state  arrangements  for  ex¬ 
tended  intervals  of  the  field,  while  A^  and  A2B  are  groundstates 
only  just  at  the  critical  field  Hc1/|j|*4,  and  a4B  as  well  as  A5B 
are  groundstates  only  right  at  the  critical  field  H  ,/|j|*12.  On 

C4 

the  other  hand,  just  at  the  critical  field  the  system  is  so  de¬ 
generate  that  it  has  a  finite  entropy.  Using  boundary  conditions 
which  do  not  favor  a  particular  ordering  leads  to  a  disordered 
system  / 1 2 , 1 7/ .  This  fact  can  be  most  easily  seen  from  the  following 
simple  argument  /12/.  Consider  two  neighboring  cells  of  the  A-^B 
Structure  (Fig.  9a)  ,  and  focus  attention  upon  the  spin  in  the  center 
of  the  plane  joining  the  cells.  Four  of  its  bonds  to  nearest  neigh¬ 
bors  are  energetically  favorable  (drawn  dash-dotted)  while  the  other  eight  are 
unfavorable  (drawn  broken).  Therefore  at  the  critical  field  H  , 

Cl 

this  spin  can  be  overturned  with  no  energy  cost.  Since  there  are 
no  bonds  between  spins  on  the  same  sublattice,  we  have  two  degenerate 
states  for  ever*  spin  on  this  sublattice  if  we  keep  the  spin  con¬ 
figuration  of  the  three  other  sublattices  fixed.  As  a  result  we 
would  have  one  sublattice  of  down  spins,  two  sublattices  of  up  spins 
and  one  sublattice  disordered,  i.e.  a  ground-state  entropy  of 
S0»Un2)/4  and  magnetization  m0«1/4  (cB-3/8)  at  the  critical  field. 

Of  course,  this  picture  of  perfect  order  on  three  sublattices  and 
perfect  disorder  on  the  fourth  is  too  simple  (yielding  however  a 
lower  bound  on  SQ).  What  we  rather  expect  to 
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happen  is  a  disordering  of  all  sublattices,  where  short-range 
correlations  are  created.  Consider  as  an  example  a  group  of 
eight  cells  (Fig.  9b)  where  four  such  loose  spins,  which  are 
nearest  neighbors  on  their  sublattice,  have  been  turned  around: 
now  it  is  a  spin  of  the  other  sublattice  in  the  "antiferro¬ 
magnetic"  plane  which  is  loose  and  can  be  turned  around.  A 
similar  degeneracy  exists  for  spins  in  the  "ferromagnetic" 
planes  of  the  A^B  structure  (Fig.  9c) .  It  is  not  clear  to  what 
extent  the  above  simple  estimates  of  SQ  and  mQ  are  modified  by 
this  more  complicated  disorder.  The  Monte-Carlo  results  at  T+0 
for  H=Hc1  are  extremely  close  to  m=1/4,  however;  thus  the 
extrapolation  of  the  phase  boundaries  between  the  disordered 
phase  and  the  two-phase  regime  ending  at  the  point  T»0,  cB=3/8 
may  have  some  significance.  We  did  not  attempt  to  extrapolate 
the  phase  boundaries  between  the  ordered  phases  and  the  two- 
phase  regimes,  however.  Studying  the  low-temperature  behavior 
by  Monte  Carlo  methods  is  rather  delicate,  as  the  convergence 
may  be  strongly  slowed  down  for  T  -*  0,  and  metastable  states 
can  hardly  be  distinguished  from  stable  ones. 

Fig.  10  presents  an  example  for  the  results  on  order 
parameters  obtained  in  the  present  study.  While  the  Cowley 
short  range  order  parameter  -a1  increases  monotonically  with 
concentration  cg  up  to  cB*1/2  at  high  temperatures, 

it  develops  a  minimum  at  lower  temperatures 

around  cB*>3/8.  This  behavior  also  reflects  the  fact  that  in 
the  vicinity  of  this  concentration  there  is  much  less  ordering 
tendency,  as  noted  above  from  ground-state  considerations.  At 
lower  temperatures  and  stoichiometric  compositions  a,  soon 
reaches  its  saturation  values  ( A,B  structure:  a , =- 1 / 3  ,  AE  struc- 
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ture:  a^-1/3),  while  it  varies  linearly  with  cg  in  the  two- 
phase  regions.  At  the  stoichiometric  compositions,  the  tem¬ 
perature  variation  of  the  order  parameters  is  found  to  be 
neither  very  sensitive  to  the  model  nor  the  approximation, 
as  Fig.  10  shows,  and  even  experimental  data  for  the  Au^Cu- 
system  / 3 3/  are  in  rough  agreement  with  the  calculation.  How¬ 
ever,  from  the  phase  diagram  of  the  present  model  (Fig.  5) 
it  is  already  clear  that  it  cannot  represent  the 
experimental  Au^Cu  system.  Thus  we  suggest  that  one  can  check 
both  the  accuracy  of  a  model  as  well  as  the  accuracy  of 
approximations  when  one  studies  off-stoichiometric  compositions. 

Fig.  11  shows  the  concentration  dependence  of  ordering 
energy  and  long  range  order  parameter.  It  is  seen  that  the 
ordering  energy  changes  very  little  throughout  the  ordered 
regime  and  even  above  the  maximum  transition  temperature  it 
is  close  to  the  variation  in  the  ground  state.  This  is  again 
a  reflection  of  the  fact  that  both  energy  differences  and 
entropy  differences  between  the  various  ordered  phases  are 
very  small,  and  hence  the  entropy  of  the  phases  needs  to  be 
calculated  very  accurately,  if  one  wants  to  make  reliable 
predictions  about  the  ordering  behavior  of  this  model  system. 

This  is  probably  the  reason  why  many  approximate  methods  give 
such  unreliable  results. 

At  low  temperatures,  the  long  range  order  parameter  of  the 
A^B  structure  has  its  maximum  at  the  stoichiometric  composition, 
but  at  higher  temperatures  it  occurs  at  cg  =  0.26  rather 
than  at  cg=0.25,  consistent  with  our  finding  that  the  maximum 
transition  temperature  cf  this  phase  occurs  at  cg*0.265  !Fig.  7). 
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We  find  that  the  transition  to  the  disordered  phase  occurs 
typically  when  the  long  range  order  parameter  has  decreased 
to  roughly  0.6  to  0.7.  This  may  perhaps  be  understood 
as  coming  from  " f rust  ration” .  The  system  can  not  maintain 
itself  in  a  state  of  long  range  order  -  there  are  too  many 
opportunities  for  spins  to  flip.  This  results  in  a  complete 
break-down  of  the  ordered  state  and  thus  in  a  first  order 
transition . 


s 
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4.  Results  for  a  Next-Nearest  Neighbor  Allov  (a=-0.25) 


Recently  Sanchez  and  de  Fontaine  /4/  applied  the  CV  method 
to  obtain  the  temperature-composition  phase  diagram  of  the  alloy 
at  a=-0.25.  This  was  our  motivation  to  perform  canonical  ensemble  Monte 
Carlo  calculations  at  c  =0.25  and  c  =0.50.  We  find  the  transitions 

Q  O 

to  be  of  first  order,  in  agreement  with  the  CV  prediction,  and 
estimate  k„T  /|j|=1.116  +  0.002  (c  =0.25),  while  the  CV  result 

D  C  D 

is  kQTc/ | J | =1 .17  /4/.  a  similar  deviation  of  only  a  few 
percent  occurs  for  co=0.5  where  we  get  k_T  /|j| =1.431  +  0.004, 
while  the  CV  result  is  /4/  k  T  /|j|=1.50.  Thus  again  the  CV 

D  C 

method  is  quite  reliable  for  the  stoichiometric  compositions. 
Preliminary  results  at  cB=0.33  indicate  that  the  CV 
method  is  less  reliable  at  non-stoichiometric  compositions  also 
in  this  case. 

Since  a  fairly  complicated  phase  diagram  has  been 
predicted  (Fig.  12),  a  more  complete  investigation  of  this 
system  by  Monte  Carlo  methods  would  be  desirable,  but  this  is 
beyond  the  scope  of  the  present  investigation. 
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5.  Summary 


We  have  given  a  detailed  description  of  the  Monte  Carlo 
method  for  the  study  of  alloy  phase  diagrams,  order  parameters 
etc. ,  and  described  results  for  the  face-centered  cubic  lattice 
with  nearest  neighbor  interaction  in  detail.  We  have  found 
that  the  phase  diagram  obtained  from  approximate  analytic 
methods  are  generally  not  very  reliable.  Only  the  cluster 
variation  method  with  tetrahedral  clusters  has  yielded  both 
transition  temperatures  and  order  parameters  reliably  (i.e.  , 
within  an  error  of  at  most  a  few  percent),  for  stoichiometric 
compositions.  Other  approximate  methods  used  so  far  are  quite 
inaccurate  even  then.  The  present  results  are  therefore  a 
challenge  to  develop  better  theoretical  methods  for  the 
calculations  of  phase  diagram  of  these  " f rustrated-type" 
systems.  Further  work  on  the  cluster  variation  method  with 
higher  order  clusters  also  seems  very  desirable.  Needless  to 
say  we  have  only  considered  a  very  simple  model,  which  is  not 
expected  to  represent  faithfully  any  real  alloy  system  -  but 
more  complicated  realistic  interactions  are  not  likely  to 
make  the  currently  available  approximation  methods  any  more 
reliable.  The  methods  are  insufficient  in  our  case  because 
the  treatment  of  configurational  entropy  of  the  alloy  is  too 
crude,  and  this  fact  will  likely  also  be  true  in  the  more 
complicated  cases.  Computer  simulation  methods  do  appear  to 
yield  numerically  accurate  information  on  both  short  and  long 
range  order  parameters  and  the  phase  diagram,  but  they  are 
expensive  to  use  because  of  the  large  amounts  of  computer  time 
required;  and  the  computations  have  to  be  repeated  for  each  set 
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of  parameters.  Nevertheless,  in  future  work  we  hope  to  apply 
the  same  methods  to  more  complicated  models  whicn  are  closer 
to  real  systems.  As  a  first  step  in  this  direction  we  have 
given  a  few  results  for  an  alloy  with  a  next-nearest  neigh¬ 
bor  interaction. 

Some  conclusions  which  emerge  from  this  study  may  stimulate 
more  experimental  work  on  ordering  alloys.  We  have  found  that 
the  behavior  of  stoichiometric  alloys  is  rather  insensitive 
to  the  parameters  of  the  model,  and  to  the  accuracy  of  the 
approximations  used.  On  the  other  hand,  observing  concentra¬ 
tion  dependences  of  long-  and  short  range  order  parameters 
(Figs.  10,12)  would  be  much  more  sensitive  tool  for  checking 
whether  a  model  faithfully  represents  a  real  system. 

There  have  also  been  developed  recently  some  rigorous 
methods  for  describing  the  phase  diagram  of  complex  alloy 
systems  at  low  temperatures.^**'^0  Unfortunately  these  methods 
are  restricted  at  the  moment  to  very  low  temperatures  and 
require  that  the  degeneracy  of  the  ground  state  be  finite. 

As  we  have  noted  this  is  not  the  case  at  the  "interesting" 
values  of  the  field  and  of  the  composition  at  least  for  the 
model  studied  here.  Fortunately  this  restriction  may  soon  be 
removed. ^ 
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Fiaure  Captions 


Fig.  1 :  a)  One  cube  of  the  fee  lattice  shown  for  the  Afl 

structure  of  CuAuI-type,  A-atoms  being  represented 
as  "spin  u;."  ,  B-atoms  as  "spin  down";  b)  two  cubes, 
A2B2  structure;  c)  one  cube,  A^B  structure  of  Cu^Au- 
type;  d)  two  cubes,  A^B  structure  of  Al^Ti-type; 
e)  six  cubes,  A2B  structure  of  NijV-type.  Numbers  in 
case  a)  indicate  the  sublattice  labelling. 

Fig.  2:  Internal  energy  (upper  part),  order  parameters 

(middle  part)  and  magnetization  (lower  part)  of 
the  nearest-neighbor  Ising  anti ferromagnet  plotted 
vs.  field  at  k  T/|j|=1.5.  Both  data  from  canonical 

O 

(C)  and  grand  canonical  (GC)  simulations  are 
included.  Points  with  arrows  denote  states  with 
unstable  ordering  relaxing  towards  the  disordered 
phase.  The  estimates  for  the  three  critical  fields 

and  H  ..  are  also  indicated, 
cl  cl  cz 

Fig.  3;  Internal  energy  (upper  part),  order  parameters  (middle 

part)  and  magnetization  (lower  part)  of  the  nearest- 
neighbor  Ising  anti ferromagnet  plotted  vs.  field  at 
kBT/|j|-1.7. 

Fig.  4:  Internal  energy  (upper  part),  order  parameter  m.  „ 

A3B 

(middle  part)  and  magnetization  (lower  part)  of 
the  nearest-neighbor  Ising  antiferromagnet  plotted 
vs.  field  at  kBT/|j|*2.0. 


Fig.  5:  Phase  diagram  of  the  nearest-neighbor  face-centered 

cubic  lattice  in  the  temperature-composition  plane 

(upper  part)  and  in  the  temperature-field  plane. 

Ordered  structures  are  indicated.  All  transitions 

are  found  to  be  of  first  order.  Open  circles  represent 

the  estimates  for  the  critical  fields  H' .  ,H".,  and  H  - 

ci  ci  cz 

mentioned  in  the  text. 


Fig.  6: 


Fig.  7: 


Fig.  8: 


Fig.  9: 


Fig.  10 


Fig.  11 


Temperature  variation  of  the  magnetization  (upper 
part) ,  order  parameter  m.  _  (middle  part)  and 

A^  o 

internal  energy  (lower  part)  for  various  values 
of  the  field.  Broken  curves  indicate  metastable 
states. 

Temperature-composition  phase  diagram  of  the  fee 
binary  alloy  with  nearest-neighbor  interaction. 

Full  curves  are  the  prediction  of  the  cluster 
variation  method  in  the  tetrahedron  approximation 
/ 22/.  Note  that  in  the  pairwise  interaction  model 
all  phase  diagrams  are  symmetric  around  cb*1/2. 

Variation  of  the  ordering  energy  at  T=0  with 
composition  (a)  /2,3/  or  field  (b) .  The  various 
ordered  structures  (or  mixtures  of  them)  are 
indicated. 

Configurations  of  neighboring  cells  in  the  A^B 
structure  where  spins  become  "loose"  at  the  critical 
field  Hci .  For  explanations  cf.  text. 

Variation  of  the  Cowley  short  range  order  (SRO) 
parameter  a1  with  composition  (upperpart)  and  plot 
of  a1  and  long  range  order  parameter  (LRO)  4/  vs. 
temperature  at  c0=1/4  (lower  part) .  Full  curves  are 
the  MC  results,  broken  curves  the  CV  method  (Ref. 

20) ,  dash-dotted  curve  the  Kittler-Falicov  theory 
(Ref.  9),  open  circles  are  experimental  data  for 
AuCUj  (Ref.  33) . 

a)  Ordering  energy  of  the  fee  nearest  neighbor  alloy 
plotted  vs.  composition  at  two  temperatures,  b)  Long 
range  order  parameter  plotted  vs.  composition  at  three 
temperatures.  Dots  denote  Monte  Carlo  results  obtained 
from  GC  calculations.  Broken  straight  lines  indicate 
the  two-phase  regime. 


i 


Fig.  12:  Phase  diagram  of  the  fee  next  nearest  neighbor 

alloy  (a=-0.25)  according  to  the  cluster  variation 
method  /4/.  Points  denote  present  Monte  Carlo 
results  for  the  transition  temperatures  at  Cg=1/4 
and  c  =1/2. 
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